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^^ ■ Abstract 

^T) ■ To implicitly restart the second-order Arnoldi (SOAR) method proposed by Bai and 

CN I Su for the quadratic eigenvalue problem (QEP), it appears that the SOAR procedure 

must be replaced by a modified SOAR (MSOAR) one. However, implicit restarts fails 
to work provided that deflation takes place in the MSOAR procedure. In this paper, we 
first propose a Refined MSOAR (abbreviated as RSOAR) method that is based on the 
refined projection principle. We derive upper bounds for residual norms of the approx- 
,J^ ■ imate eigenpairs obtained by the MSOAR and RSOAR methods. Based on them, we 

j^ I propose a reliable tolerance criterion for numerical breakdown that makes the MSOAR 

and RSOAR methods converge to a prescribed accuracy. This criterion also serves to 
decide numerical deflation. We consider the central issue of selecting the shifts involved 
when implicitly restarting the MSOAR and RSOAR algorithms. We propose the exact 
and refined shifts for the two algorithms, respectively, and present an effective approach to 
treat the deflation issue in implicit restarts, so that the implicit restarting scheme works 
^-j- I unconditionally. Numerical examples illustrate the efficiency of the restarted algorithms 

OS ■ and the superiority of the restarted RSOAR to the restarted MSOAR. 
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1 Introduction 

Consider the large quadratic eigenvalue problem (QEP) 

{X^M + XC + K)x = 0, (1) 

where M, C, K are nxn real or complex matrices with M nonsingular. QEP's arise in many 
applications, see, e.g., [21 [27]. We are interested in a few largest eigenvalues in magnitude or 
a few eigenvalues nearest to a target in the complex plane. A commonly used approach is 
to linearize the QEP. There are a number of linearizations available |27j . For example, if M 
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is invertible, one of the most often used linearizations is to transform it to the generahzed 
eigenvalue problem 
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which can be further reduced to the standard linear eigenvalue problem 
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where A = —M ^C,B = —M ^K. Note that ([3]) corresponds to the monic QEP 

(A^ -XA- B)x = 0. (4) 

The mathematical theory on ([2]) and ([3]) has been well established and a number of nu- 
merical methods have been available for solving it [H [6l [23l [25l [28] . A serious drawback of 
the linearization is that general numerical methods do not take the structures of ([2]) and 
([3]) into account, so that resulting projected matrices do not possess the same structures as 
those of the original problems. Moreover, essential structures of the QEP, such as possible 
symmetries of M, C and/or K are not preserved. For example, if M, C and K are real 
and M is symmetric positive definite, K is symmetric and C is skew-symmetric, then the 
eigenvalues of the QEP appear in pairs (A, —A) if A is real or purely imaginary and in quadru- 
plets (A, —A, A, —A) if A is complex. However, the eigenvalues of projected matrices for the 
linearized problems do not have such spectrum structures. Therefore, it is essential to work 
on the QEP directly. We will develop methods that can not only preserve essential structures 
and properties of the QEP but also take some advantages of the Arnoldi type methods for 
the linearized problem. 

The Second-Order Arnoldi (SOAR) method proposed by Bai and Su [3] is a Rayleigh-Ritz 
method that works on ([1]) directly. They propose a SOAR procedure that first computes an 
orthonormal basis for a second-order Krylov subspace generated by the matrices A and B 
simultaneously in an elegant way and then projects ([T]) onto this subspace to get a projected 
QEP that preserves the physical structure of the original QEP. They have established some 
relationships between the SOAR method and the Arnoldi method for the linear problem ([3]) 
and shown that it has some merits of the Arnoldi method. The SOAR procedure is also used 
in dimension reduction of large scale second-order systems [4] . It preserves the second-order 
structure and achieves the same order of approximation as the standard Arnoldi method 
via linearization. A unified general convergence theory has recently been established in 
[9], which has generalized some of the results on the Rayleigh-Ritz method for the linear 
eigenvalue problem. It is proved that for a sequence of projection subspaces containing 
increasingly accurate approximations to the desired eigenvector there is a Ritz value that 
converges to the desired eigenvalue while the Ritz vectors converge conditionally and may 
fail to converge. Alternatively, we can compute a refined Ritz vector whose unconditional 
convergence is guaranteed. 

Due to the storage requirement and/or computational cost, to be practical, it is generally 
necessary to restart the SOAR method. Unfortunately, it appears that the implicit restart- 
ing technique [24] is not applicable to the SOAR method. This is a major disadvantage. 
Meerbergen [21] proposes a Quadratic Arnoldi (Q-Arnoldi) method that is an alternative 
of the SOAR method. The Q-Arnoldi method exploits the structure of the linear problem 
to reduce the memory requirements by about a half and can compute a partial Schur form 
of the underlying linearization with respect to the structure of the Schur vectors, so that 
the implicit restarting technique can be applied. Otto [22] considers some theoretical and 



numerical aspects of the SOAR method. Since implicit restarting cannot be applied to the 
SOAR method, he proposes a modified second-order Arnoldi (MSOAR) procedure that re- 
places the original special starting vector by a general one. Based on the MSOAR procedure 
and under the restrictive assumption that there is no deflation in the MSOAR procedure, 
he has developed an implicitly restarted MSOAR algorithm. It is known from [31 [22] that 
deflation or breakdown may take place in the SOAR and MSOAR procedures but deflation 
has a completely different consequence from that of breakdown. Bai and Su [3] proves that 
the SOAR method will find some exact eigenpairs of QEP ([!]) if breakdown occurs but no 
eigenpair is found if deflation takes place. A remedy strategy is given in [3] that can continue 
the SOAR procedure until the SOAR method converges. The strategy has been adapted to 
the MSOAR procedure [22]. However, it is not clear what will happen when deflation occurs 
in the MSOAR procedure. Relationships are also not known between the modified second- 
order Krylov subspace and its related standard Krylov subspace and between the MSOAR 
procedure and the standard Arnoldi process. As in the SOAR method, it turns out that 
similar relationships are basic for understanding the MSOAR methods. 

There are explicit and elegant expressions on the residual norms of approximate eigenpairs 
and the Arnoldi process, which are used to design reasonable and reliable criteria for numerical 
breakdown [231 l25l I28[ . However, no results have been established for residual norms of 
the approximate eigenpairs. Such kind of explicit expressions are certainly appealing and 
extremely important to better understand and implement the SOAR and MSOAR methods. 
Without them, one would not be able to propose reliable and reasonable criteria for deciding 
numerical breakdown. Such kind of results also appear to play a vital role in designing a 
reasonable tolerance criterion for numerical breakdown and deflation. 

There are a number of unsolved important problems in the implicitly restarted MSOAR 
algorithm. As mentioned above, the implicitly restarted algorithms proposed by Otto [22] 
critically requires that no deflation occurs in the MSOAR procedure. Whenever deflation 
takes place, implicit restarting fails completely and cannot work. This requirement limits 
the applicability and generality of the algorithms. Another central problem is reasonable 
selection of the shifts involved. Otto [22] proposes exact second-order shifts for use within 
the algorithm. Unlike those implicitly restarted algorithms for the linear eigenvalue problem, 
it is distinctive that candidates for shifts exceed the number of shifts allowed. This makes 
a selection of shifts more subtle and complicated than that for implicitly restarted Arnoldi 
type algorithms for linear eigenvalue problems. Otto does not show how to reasonably select 
'correct' shifts among the candidates. As it will appear, a reasonable selection of shifts is 
mathematically nontrivial, and we must be careful and do more work. 

In this paper, combining the SOAR and MSOAR procedures with the refined projection 
principle [HI [H] (see also [H [251 128] ) , we propose a refined second-order Arnoldi (RSOAR) 
method that can guarantee the unconditional convergence of refined Ritz vector when a 
subspace is accurate enough [9] . We will further investigate second-order Krylov subspaces 
and highlight on the defiation and breakdown issues in the MSOAR procedure. We establish 
some basic connections between the second-order Krylov subspaces and related standard 
Krylov subspaces, showing the equivalence between them as well as between the MSOAR 
method and the standard Arnoldi method. We prove that some exact eigenpairs are found 
when breakdown occurs in the MSOAR procedure. We also prove that when deflation takes 
place the first time at some step in the MSOAR procedure, the standard Arnoldi process 
must do not break down. We the focus on criteria on numerical deflation and breakdown. 
Such criteria are extremely important and serve for two purposes: one of them is to ensure 
the convergence of the MSOAR and RSOAR methods when numerical breakdown occurs; 
the other is for the numerical stability of the MSOAR procedure. To this end, we derive 
upper bounds for residual norms of the approximate eigenpairs obtained by the MSOAR and 



RSOAR methods and establish some important relationships between them and numerical 
breakdown. Based on these results, we propose a reasonable and reliable tolerance criterion 
for deciding numerical deflation and breakdown. 

After the above is carried out, we consider two key issues when implicitly restarting the 
MSOAR and RSOAR algorithms: selection of shifts and treatment of numerical deflations. 
Just as the mechanism for those implicitly restarted Krylov subspace algorithms for linear 
eigenvalues problems and SVD problems [HI [HI [T71 [181 [21], it turns out that a proper 
selection of shifts involved is one of the keys for success of the implicitly restarted algorithms. 
Based on the results on the equivalence between the second-order Krylov subspaces and 
the standard Krylov subspace as well as the MSOAR procedure and the standard Arnoldi 
process, we will see how shifts should be chosen. We propose new exact shifts that are 
different from Otto's and the reflned shifts for respective use within the implicitly restarted 
MSOAR and RSOAR algorithms. The refined shifts are based on the reflned Ritz vectors 
and theoretically better than the exact shifts. As was pointed out previously, now both exact 
and reflned shift candidates are more than the shifts allowed. This makes a selection of 
shifts subtle and complicated. We show how to reasonably select the desired shifts among 
them for each algorithm. We present an efficient algorithm to compute the exact and refined 
shift candidates reliably. Finally, we consider the critical issue of how to realize implicit 
restarts when encountering numerical deflations. We propose an effective approach to cure 
numerical deflations, so that implicitly restarted MSOAR and RSOAR algorithms are of 
general applicability and can be run unconditionally. 

The rest of this paper is organized as follows. In Section [2] we review second-order Krylov 
subspaces and the SOAR and MSOAR procedures and describe the SOAR and MSOAR 
methods. We establish a number of basic and important properties that connect second-order 
Krylov subspaces to standard Krylov subspaces and the MSOAR procedures to the standard 
Arnoldi process. In Section [3l based on the SOAR and MSOAR procedures, we propose the 
RSOAR method and discuss it in some detail. In Section [U we establish upper bounds for 
the relative residual norms of approximating (reflned) Ritz pairs. With them, we propose 
a reasonable criterion for deciding numerical deflation and breakdown. We then develop 
implicitly restarted RSOAR and MSOAR algorithms and discuss how to select reasonable 
shifts for each algorithm and propose the exact shifts and the reflned shifts for the MSOAR 
and RSOAR algorithms, respectively. We present an effective approach to treat numerical 
deflations. In Section [5l we report numerical experiments to illustrate the efficiency of the 
restarted algorithms and the superiority of the refined algorithm. 

Throughout the paper, we denote by || • || the spectral norm of a matrix and the 2-norm of 
a vector, by || • ||i the 1-norm of a matrix, by / the identity matrix with the order clear from 
the context, by the superscripts T and * the transpose and conjugate transpose of a matrix, 
by C the complex space of dimension k and by T^C^+i)^^^^ the set of {k + 1) x k real matrices. 
We denote by o"min(-^) the smallest singular value of a matrix F and by the Matlab notation 
A{i : j,k : I) the submatrix consisting of columns i to j and rows /c to / of A. 

2 Second-order Krylov subspaces, the SOAR procedure, the 
MSOAR procedure and the SOAR and MSOAR methods 

Bai and Su [3] introduce the following concepts. 

Definition 1. Let A, B he matrices of order n and the vector u ^ Q, and define 

ro = u, 

ri = Aro, 

rj = Arj_i + Brj^2 for j > 2. 



Then ro,ri,r2, . . . ,rfc_i is called a second-order Krylov sequence based on A,B and u and 
Gk{A, B; u) = span{ro, n, r2, . . . , ri^_i} a kth second-order Krylov subspace. 



Note that ([3]) is a linearization of (JH). Define the matrix 
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(5) 



of order 2n. For a 2n dimensional starting vector v, we can generate a Krylov subspace 
JCk{H,v) = span{v,Hv,H'^v, . . . ,H^~^v}. Particularly, if we choose v = [u^ , 0]^, we have 



r_i-i 



Wv, j >0 with r_i = 0. 



(6) 



This indicates that the upper half part of ICk{H, v) is just Gki^i B; u) and its lower half part 
is just Qk-i{A,B;u), respectively. These fundamental facts can be compactly expressed as 



gl_,{A,B;u) C ICkiH,v) C gl{A,B;u), 



(7) 



where g^{A,B;u) is the subspace generated by 
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see also Theorem 2.4 in [22]. Due to the equivalence of QEP (j3|) and the eigenproblem of 
H, if the eigenvector [Xx ,x ] is contained in }Ck{H,[u ,0] ), then the eigenvector x of 
QEP (JH) is contained in ^fc(^,i?;ii). Therefore, continuity arguments suggest that if there 
is a good approximation on x contained in gk{A,B;u) then there is an approximation with 
comparable accuracy on [Xx'^,x'^]'^ in }Ck{H, [m-^,0]"'"). These illustrate that for QEP (JH the 
subspace Qk{A,B;u) provides the same information as lCk{H,v) for the eigenvalue problem 
of H. This motivates us to directly solve QEP ([1]) based on Qk{A, B; u) rather than solve the 
eigenproblem of H based on ICk{H,v). 

Bai and Su [3] propose the following procedure for computing an orthonormal basis 
{Qj}'j=i of Gk{A, B;u) and an auxiliary vector sequence {pj}- 

Algorithm 1. SOAR procedure 
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Qi = u/\\u\\ 
Pi = 

for j = 1,2, ... ,/c do 
r = Aqj + Bpj 
s = qj 

for 1 = 1,2,... ,j do 
tij = q*r 
r = r- qitij 

S — S PiXij 

end for 

tj+ij = M 

if tj+ij = 0, stop 

Qj+i = r/tj+ij 
Pj+i = s/tj+ij 
end for 
This algorithm leads to the following basic results [3]. 



Theorem 1. Define Qk = [qi,q2,---, Qk] and Pk = [pi,P2, ■■■,Pk] and fk=, ^ 

I tk+lkS-k 

span{Qk} = gk{A,B-u) (8) 

and the SOAR decomposition 
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where Qk+i = [Qk,qk+i\, Pk+i = [Pk,Pk+i]- 

Q is an Arnoldi like decomposition on H starting with v = [qi,0]'^, but the columns of 

Qk 

Pk 

are a non-orthogonal basis of ICk{H, v). Furthermore, qi,q2,---,qk form an orthonormal basis 
of Qk{A, B; u) and p2, ■ ■ ■ ,Pk form a non-orthogonal basis of the subspace Gk-ii^y B] u). 

If Algorithm [1] stops prematurely at step j, there are two possibilities. One possibility 
is that Tj, i = 0, . . . ,j are linearly dependent but [rf , r'[_-^]'-'" , i = 0,. . . ,j with r_i = are 
not. In this case, Qj^i{A,B;u) = Qj{A,B;u) but ICj+i{H,v) ^ ICj{H,v), so the Arnoldi 
process on H does not terminate at step j; see Lemma 2.4 of [22]. This situation is called 
deflation. The other possibility is that both vector sequences {rj} and {[r'[,rf_i]'^} are 
linearly dependent at step j. This situation is called breakdown and the Arnoldi process on 
H terminates; see Theorem 3.1 of [3]. If deflation occurs in the SOAR procedure at step 
j, JCj{H,v) does not contain any eigenvector of H, which implies that gj{A, B;u) does not 
contain any eigenvector of QEP (j4]); if breakdown occurs in the SOAR procedure at step j, 
Qj{A, B; u) contains 2j exact eigenvectors of QEP @. Therefore, deflation must be remedied. 

Bai and Su [3] present the following algorithm that can remedy deflation. 

Algorithm 2. SOAR procedure with deflation 
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qi = 'u/ll^ll 

for j = 1,2,... ,k do 

r = Aqj + Bpj 
s = qj 

for 1 = 1,2,... ,j 
tin = q-r 

^ qitij 
Pitij 



do 



r = 

s = s - 
end for 

tj+ij = IMI 
if tj+ij = 

if s £ span{pi\7 

break 
else deflation 
reset tj+ij = 
qj+i = 

Pj+i = s 
end if 
else 

Qj+i = r/tj+ij 

Pj + l = s/tj + ij 

end if 
end for 



(Zi = 0, 1 < i < j} 



In the above procedure, if ij+ij = at line 12 and deflation occurs, we simply take 
Qj+i = and set tj+ij to one. To decide if s € span{pi\i : qi = 0,1 < i < j}, the Gram- 
Schmidt orthogonalization with iterative refinement can be used, as suggested in [31 E2]- The 
procedure can continue until it breaks down. When deflation occurs, the nonzero vectors in 
the sequence {qi} still span the second-order Krylov subspace Qk{A, B;u) whose dimension 
is smaller than k. For details on Algorithm [2l we refer the reader to Bai and Su [3], where 
they present a SOAR procedure with deflation and memory savings by exploiting pi = 0. 

For the SOAR procedure with deflation, Theorem[T]and relation Q are still true. Because 
of ([7]), we can easily justify the following result. 

Theorem 2. If the Arnoldi process on H breaks down at step k, Algorithmic breaks down at 
step k — 1; if Algorithmic breaks down at step k, the Arnoldi process breaks down at step k. 



As was realized |21ll22j . a serious disadvantage of the SOAR procedure is that the implicit 
restarting technique cannot applied to it since the updated pi is not zero any more when 
implicitly restarting it. Otto |22] proposed modifying the original starting vector pi = 
in Algorithms [IHl] and allowing it to be nonzero. This leads to a modified second-order 
Krylov subspace and a modified second-order Arnoldi (MSOAR) procedure that generates 
an orthonormal basis of it, as shown below. 

Definition 2. |22j Let A and B be n x n matrices and ui,U2 G C" nonzero vectors, and 
define the sequence 

ro = ui, 

ri = ArQ + Bu2, 

rj = Arj^i + Brj^2 for j > 2. 

Then ro, ri, r2, . . . , rk-i,U2 is called the kth second-order Krylov sequence based on A, B and 
ui,U2, and QkiA, B;ui,U2) = span{rQ,ri,r2, ■ ■ ■ ,rk-i,U2} the kth second-order Krylov sub- 
space. 

Note that the subspace QkiA, B; ui,U2) is spanned by the upper block of the sequence 
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n = Aro + Bu2 



rk-i = Ark~2 + -Brfc_3 
rk-2 



U2 





(10) 



It is clear that the subspace spanned by its upper block vectors contains the subspace spanned 
by its lower block vectors. For }Ck{H, v) with v = [uj, u'^]'^, we easily see its upper and lower 
half parts are contained in Qk{A, B; ui,U2) and Qk-i{A, B; ui,U2), respectively, and vice versa. 
Similar to ([7]), these simple but fundamental facts can be expressed as 



gl^^{A,B;ui,U2) C ICkiH,v) C gliA,B;ui,U2). 



(11) 



Therefore, for QEP (j4]), Qk{A,B;ui,U2) provides the same information as ICk{H,v) for the 
eigenproblem of H. 

Algorithm [3] describes the MSOAR procedure [22] that can remedy defiation and gener- 
ates an orthonormal basis of QkiA, B; ui, U2) of dimension k -\- 1. 

Algorithm 3. MSOAR procedure 

Input: Matrices A, B, nonzero vectors ui, U2, and steps k. 

Output: // no premature stop occurs, the nonzero vectors among qi, . . . ,qk+i form, an or- 
thonormal basis for Qi:{A,B;ui,U2), the sequence {pi, ■ ■ ■ ,Pk} spans ty^_i(A, i?; ni,n2) and 



A; + 1) X k upper Hessenberg matrix T^ 

'q^ = ^,pi = ^, / = 0. 
||«i|| ' ^^ ll"i|| 

for j = 1,2,..., k do 
r = Aqj + Bpj 



[tij] is returned. 
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for i=l,2,...,j do 
tij = q*r 
r = r- tijQi 

S — S t^ijPi 

end for 

tj+ij = Ikll 
if tj+ij = 

1 = 1 + 1, tj+ij = 1, qj+i = 0, pj+i = s, fi= pj+i. 

if (/ = Okpj+i G span{qi, . . . , qj,pi}) or {I > O&ipj+i G span{fi, ..., fi-i}) 
break, go to step 21. 

end if 
else 

Qj+i = r/tj+ij 

Pj+i = s/tj+ij 
end if 

end for 

qj+2 = Pi 

for i = 1 : j + 1 do 

7 = qiqj+2 
qj+2 = qj+2 - iqi 

end for 

qj+2 

end if 



gj+2 

Iki+all 



Like the SOAR procedure, the MSOAR procedure may stop prematurely at step j < k. If 
so, deflation or breakdown will take place. That is, if tq, ri, . . . , Tj_i , U2 are linearly dependent 
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are linearly independent, then deflation occurs; if 



both vector sequences are linearly dependent, then breakdown occurs. 

From Algorithm [3] and the key relation (jlip . we can easily establish the following results, 
which are analogues of Theorems [TH2l 

Theorem 3. The columns of Qk+i generated by the k-step MSOAR procedure form an 

Qk 



orthonormal basis of Qii.{A,B;ui,U2), ([9]) holds and the columns of 



Pk 



form a non- 



orthogonal basis of lCk{H,v). Furthermore, if the Arnoldi process on H breaks down at step 
k, then Algorithmic breaks down at step k — 1; if Algorithmic breaks down at step k, the 
Arnoldi process on H breaks down at step k. 

This equivalence theorem indicates that breakdown is a lucky event since there are 2k 
exact eigenvectors of QEP ([1]) contained in Qk{A,B\ui,U2) if Algorithm [3] breaks down at 
step k and no deflation occurs before step k. However, if deflation occurs, the situation 
becomes completely different, as the following results reveal. 

Theorem 4. Assume that I is the smallest integer such that 

gi+i{A,B;ui,U2) = Gi{A,B;ui,U2). 

Then we must have 

JCi+i{H,v)^ICi{H,v). 



Proof. From the definition of ICi{H,v), tlie assertion amounts to proving tliat 
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(12) 



holds if and only if 70 = 71 = • • • = 7; = 0. Under the assumption on Z, we know vq, . . . ,ri,U2 
are linearly dependent but ro, . . . , r;_i, ii2 are linearly independent. So 

i-i 

y ajTj + aiU2 = means that oq = ai = ■■■ = «/ = 0. (13) 

j=0 

Since (1121) can be written as 



70^0 + 71^1 H ^ 7«^/ = 0, 

7o^i2 + 71^0 H h 7in-i = 0, 

from (I13p we must have 70 = • • • = 7i = 0. Therefore, the assertion is proved. D 

An important consequence of this theorem is that the first premature stop of Algorithm [3] 
must be deflation rather than breakdown and no exact eigenvector of QEP ([T]) is contained 
in Qi{A, B; ui, U2), so we must continue the algorithm for more steps. 

The SOAR and MSOAR methods are orthogonal projection methods. They project ([!]) 
onto Qi^{A, B;u) and Gk-i{A, B;ui,U2), respectively, and solves a, k x k projected QEP 

{dH'h + eCk + Kk)g = 0, (14) 



where M^ = QlMQk, Ck = QlCQk and K^ = QlMQk- It is clear that (^ preserves 
essential structure of the original QEP ([T]), such as possible symmetries, skew-symmetries, 
positiveness of M, C and K as well as possible spectrum symmetries. Suppose that the 
{9,g) are the eigenpairs of ()14p . Then the SOAR or MSOAR method uses the Ritz pairs 
{9,y = Qkd) to approximate some of the eigenpairs of ([1]). The 9 and y are called the Ritz 
values and Ritz vectors of ([I]) with respect to Qk{A, B; u) or Qk~i{A, B; ui,U2)- If Algorithm [3] 
breaks down at step k and no deflation occurs before step k, the MSOAR method will find 
at least k exact eigenpairs of ([T]), as shown in J9]. 

3 A refined second-order Arnold! (RSOAR) method 

As is well known, for a sequence of subspaces containing increasingly accurate approximations 
to the desired eigenvectors, orthogonal projection methods may fail to compute eigenvectors 
|101 119|. To correct this deficiency, a refined projection principle is proposed in |11|, [T3] (see 
also [251 [28] ) for the linear eigenvalue problem that extracts the best approximate eigenvectors 
from a given subspace in the sense that the residuals formed with certain approximate eigen- 
values available are minimized over the subspace. The refined projection methods correct 
possible non-convergence of eigenvectors and computationally viable. The refined projection 
principle has been extended to solve the SVD problems [El [171 [181 [20] . Since the SOAR and 
MSOAR methods are orthogonal projection (Rayleigh-Ritz) methods, they also have the 
possible non-convergence for computing eigenvectors, as has already been shown [9]. To this 
end, the refined projection principle is extended to the QEP that finds a certain refined Ritz 
vector whose convergence is guaranteed when the projection subspace is accurate enough [9]. 
Below we combine the SOAR procedure with the refined projection principle and propose a 
refined SOAR or MSOAR method. 



Suppose that we have computed an approximate eigenvalue 9 by the SOAR or MSOAR 
method. The Refined SOAR or MSOAR method, written as the RSOAR method hereafter, 
seeks a unit length vector u € Qk{A^B]u) or Qk_i[A,B]Ui,U2) satisfying the optimal re- 
quirement 

\\{9^M + eC + K)u\\= min \\{e'^M + OC + K)u\\ (15) 

u e gk{A,B;u) or gk-i{A,B;ux,U2) 
\\u\\ = 1 

and uses it as an approximate eigenvector, called the refined Ritz vector or more generally 
a refined eigenvector approximation. Since the columns of Q^ form an orthonormal basis of 
Qk{A,B;u) or Qk-i{A,B]Ui,U2), (fTSJl amounts to seeking a unit length vector z € C'^ such 
that u = QkZ with z the solution of 

min \\{e'^M + eC + K)Qkz\\. (16) 

Obviously, z is the right singular vector of the matrix O^MQf^ + 9CQk + KQ}^ associated 
with (Tmin(^^Af(5A: + 9CQk + KQk). However, the direct computation of the SVD may be 
expensive. Assume that the matrix is real and k <ti n. If is real, the cost of Golub-Reinsch's 
SVD algorithm is about Ank'^ flops and that of Chan's SVD algorithm is about 2n/c^ flops [6l 
p. 254]. Suppose we want to compute m eigenpairs with m < k. The total costs are Anrnk"^ 
and 2nmk'^ flops, respectively. 

The first author in [15j has proposed a cross-product matrix based algorithm for comput- 
ing the SVD of a matrix, which can be much more efficient than the standard SVD algorithms. 
Applying the algorithm to our case, we form the cross-product matrix 

Bk = {9^MQk + 9CQk + KQkY {9^MQk + 9CQk + KQk) , 

which is the symmetric (Hermitian) (semi-) positive definite and whose smallest eigenpair 
is ((T^;^(0^MQfc + 9CQk + KQk),z). We then compute the eigensystem of B^ by the QR 
algorithm to get z. In finite precision arithmetic, it is proved in [15] that the computed 
eigenvector is an approximation to z with accuracy O(emach) and the square root of the 
Rayleigh quotient of Bj. with respect to the computed eigenvector is an approximation to 
c^mm{9'^MQk + 9CQk + KQk) with accuracy O(emach) provided that the second smallest 
singular value of 9'^MQk + 9CQk + KQ^ is not close to the smallest one, where emach is the 
machine precision. 

Let us now look at the computational cost of this algorithm. Define 

Wi = MQk, W2 = CQk, Ws = KQk, 

which are available when forming the projected QEP and do not need extra cost. Then 

Bk = \9\'^ WiWi+ I 9 |2 W2W2 + 1^3*1^3 + 99^W^W2 + 99^W2Wi (17) 

+PW^W3 + 9^W2Wi + 9W2W3 + 9W^W2, 

where the bar denotes the complex conjugate of a scalar. Assume that Wi,W2 and W3 are 
real and note that Bk is Hermitian for a complex 9 and real symmetric for a real 9. Then we 
only need to form the upper (lower) triangular part of i?fc, which involves the upper (lower) 
triangular parts of the nine matrices W*Wj, i,j = 1, 2, 3. So the total flops are about 9n/c^. 
With W*Wj, i,j = 1,2,3, we only need 0{k^) flops to form Bj^ for either a real or complex 
9, negligible to 9nk'^ flops. So we only need 9n/c^ flops to form in Hermitian matrices B^s 
for m approximate eigenvalues 0's. We then compute the complete eigensystems of these m 
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-Bfc's by the QR algorithm using 0{mk^) flops. This means that we can compute m right 
singular vectors z's using about 9nk'^ flops when mk <C n, a natural requirement in practice. 
In view of the above, we see that this cross-product based algorithm is more efficient than 
Golub-Reinsch's SVD algorithm when m >3 and Chan's SVD algorithm when m > 5. 
We can now present a basic (non-restarted) RSOAR method. 

Algorithm 4. The RSOAR algorithm 

1. Given the starting vector u or ui,U2, run the SOAR or MSOAR procedure to generate 

an orthonormal basis Qk of Gk{A,B;u) or Qk~l{A,B■,Ul,U2)■ 
2. Compute Wi = MQk, W2 = CQk and W3 = KQk- 

3. Compute Mk = QlWi, Ck = QIW2 and Kk = QIW3, solve the projected QEP 

{efMk + 9iCk + Kk)gi = (18) 

and select m Ritz values 9i as approximations to the m desired eigenvalues Aj. 

4. For each chosen 6i, i = 1,2, ... ,m, form B^ and compute the eigenvector zi of Bk 
associated with its smallest eigenvalue by the cross-product matrix based algorithm. 

5. Test accuracy of the m refined eigenpairs {Oi,Ui) by computing the relative residual 
norms 

\\{9fM + eiC + K)Qk~z^\\ 

||M||i + ||C7||i + ||K||i ' ^""^ 

where \\ ■ ||i is the 1-norm. 

4 Implicitly restarted algorithms 

This section consists of four subsections. In Section [4.11 we consider how to determine defla- 
tion and breakdown numerically. We propose a reasonable criterion for deciding numerical 
deflation and breakdown. In Section 14.21 we briefly show why implicit restarting cannot be 
applied to the SOAR procedure and the MSOAR procedure is needed instead. Under the as- 
sumption that no deflation occurs, we review how to implicitly restart the MSOAR procedure 
and develop implicitly restarted MSOAR algorithm and RSOAR algorithm. In Section 14.31 
we discuss how to select shifts and propose exact and refined shifts for two algorithms, re- 
spectively. In Section 14.41 we propose an effective approach to treat the deflation issue in 
implicit restarts, so that implicitly restarted algorithms can be run unconditionally and are 
of generality. 

4.1 Determination of numerical deflation and breakdown 

We consider the crucial issue of how to decide deflation and/or breakdown of Algorithm [3] 
numerically. In line 13, we use the Gram-Schmidt orthogonalization with refinement to de- 
termine if Pj+i is in span{qi, . . . , gj, M2} or span{fi, . . . , fi-i}. In finite precision arithmetic, 
tj+ij is rarely zero exactly and, in general, can only be numerically small. If tj+ij satisfies 

*-'+'^' < tol (20) 



l^lli + ll<^lli + ll-^^"lli 



for some suitably small tolerance tol, we accept it as zero numerically. A reasonable choice for 
tol should satisfy the two requirements as stated in the introduction of Section 2) tol should 
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not be too small to cause the instability of Algorithm [3] since a very small tol may cause 
great growth in Pj+i in line 18; tol should make the MSOAR and RSOAR methods converge 
within a user-prescribed tolerance whenever numerical breakdown occurs. To meet these two 
requirements, we will prove that tol should be as small as the convergence tolerance. 

If (pUj) is satisfied, numerical deflation and/or breakdown may take place. We decide 
whether or not it is numerical deflation in the following way: In line 13 of Algorithm [3l 
before the Gram-Schmidt orthogonalization is run, we first normalize /i,- ••,/;-! to have 
unit length. Let p be the resulting vector of orthogonalizing Pj+i against qi, . . . ,qj,U2 or 
fi, . . . , fl^i. Then if ||p|| > tol, we accept it as numerical deflation; otherwise, numerical 
breakdown occurs. 

Observe that the projected QEP (fT4]) of the original large QEP ^ over span{Qk} is 
equivalent to the generalized eigenvalue problem 



" -Ck 


-Kk' 


(Oa\ 


= e 


' Mk 


■ 


( Og 


I 





V 9 J 







/ 


V 9 



which is the projected problem of 
orthonormal matrix 



over the subspace spanned by the columns of the 



Q2k 



Qk 
Qk 



where Qk is the matrix deleting zero column(s) of Qk- The projected problem is further 
equivalent to the standard linear eigenvalue problem 



-Mk'Ck 
I 



-Mk'Kk 




9 



9 



Recall that the k-step MSOAR procedure can be written as 



A B 
I 



Qk 
Pk 

Qk+2 



Qk 
Pk 

{I- 



Tk + tk+lk 

Qk+iQl+i)pi 



Qk+l 
Pk+1 



"k^ 



\\{i-Qk+iQU,)p 



with Tk = \t 
projection matrix of H 



jjj and Ck the kth coordinate vector of dimension k. Note that Tk is the oblique 

'A B 1 

with the right subspace span{[Q'^ , P^]'^} and the left 



/ 



1,2,..., 

Qk 
Pk 



k of Tk are the Petrov-Ritz 
s,- with Si the normalized 



subspace span{[Q'^,0]^}. So the eigenvalues Ui, i = 
values and the corresponding Petrov-Ritz vectors Wi 

eigenvectors of Tk associated with the eigenvalues Vi. 

Specially, if [Q'k,Pk]'^ is column independent and the MSOAR procedure breaks down 
at step k, i.e., ifc+ifebfe+iiPfc+i]"^ = 0, then span{[Q'^ , Pj[]'^} is an invariant subspace of H. 
It thus follows from Theorem [3] that the Arnoldi process on H starting with [qf ,Pi] breaks 
down at step k and the MSOAR or RSOAR method finds 2k exact eigenpairs of ([T]), where 
k is the number of nonzero columns in Qk- 

Since span{[Qjl , Pfl^]'^} C span{Q2k}, the Ritz pairs {6, [Oy^ ^y^Y) oi H over span{Q2k} 
are generally more accurate than those of H over span{[Q^ , P^f"} [lOl [19] and thus more 
accurate than the above Petrov-Ritz pairs (i^, w) [231 p. 141-2] and [281 P- 114-5]. We 
should comment that speaking of meaningful approximate eigenpairs only makes sense when 
span{[Q^, PkV^} contains good approximations to the desired eigenvectors [lOl [TH [25] . Based 
on these, we now make the following reasonable assumption. 
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Assumption: span{[Q^ , P^] } is a good subspace and the Ritz pair {9,[6y ,y ] ) is al- 
ways at least as good as Petrov-Ritz pair (u, w) when approximating the same desired eigen- 
pair (A, x), i.e., the residual norm of the former is no more than that of the latter. 



Theorem 5. Assume that {9, [9y'^,y^]^) is more accurate than {u,w) in the sense of residual 
norm, and set 

Then for the MSOAR and RSOAR methods we have 



\\i9^M + 9C + K)y\\ < Cktk+ik\e*ks\, 
\\{9^M + 9C + K)u\\ < Cktk+ik\els\. 



(22) 
(23) 



If numerical breakdown occurs at step k in the sense of (j20p . then for the MSOAR and 
RSOAR methods we have 



\\{9^M + 9C + K)y\\ 

Wh + WhTWh 

\\{9'^M + 9C + K)vi\\ 
||M||i + ||C||i + ||i^||i 



< Cfc|e^s| • tol, 

< Cfc|e^s| • tol. 



(24) 

(25) 



Proof. Prom (j2ip and A = —M ^C, B = —M ^K, the residual of (I'jw) as an approximate 
eigenpair of ([2]) is 



-C -K 
I 



w 



\w\ 



M 
/ 



w 



\w\ 



tk+lk 



M 
/ 



9fe+i 
Pk+i 



CuS 



\W\ 



(26) 



Here we normahze w by its norm \\w\\. By ||s|| = 1 and the definition oi w, we have 

\\w\\ = {\\Qksr + \\Pksff/' = {l + \\Pksff/\ 
Therefore, 

Ik.ll = tk^ik\e,s\ ^^-p^^p)T7^ • (27) 

We now investigate the residual rg of {9,{9y^ ,y'^)^) with y = Qi^g as an approximate 
eigenpair of ([2]) and seek the relationship between rg and the residual of {9, y) as an approx- 
imate eigenpair of ([1]). We have 



re 



-C -K 

I 



9y 

y 



?2m + 9C + K)y 




/(|^|2 + 1)1/2 _^ 

/{\9\'' + lf'\ 



M 



9y 

y 



/i\t 



+ 1 



,1/2 



Here we normalize {9y'^,y'^)'^) by its norm (|^p + 1)^'^. Therefore, it exactly holds that 



\\{9''MQk + 9CQk + KQ, 



\\{9^M + 9C + K)y\\ = {\9\^ + 1) 



1/2 I 



re\ 



(28) 



Since {9, [9y'^,y'^]'^) is assumed to be more accurate than {v,w) in the sense of residual 
norm, i.e., 

heW < \\ru\\, 
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we get from ([27]), (^ and i^ that 

ni /\/rii2 _i_ i[„ i[2U/2 

(1 + ||Pfcs|r)^/^ 
= Cfcifc+ifc|efcs|, 

from which and ()20p it follows that ()24p holds. 

For the RSOAR method, since the residual norm of the refined approximate eigenpair 
{6, u) is always no more than that of the MSOAR method, (I23p and ()25p are direct. D 



(j22p and ()23p are a-posteriori computable bounds for residual norms of {9,y) and {0,u). 
We should point out that the upper bounds in (j22p . (j23p . (j24p and (|25p may be conservative 



since span{[Q|^,P^]-'"} C span{Q2k} may make ||r5)|| < ||r,^|| considerably. 

Now we look into c^ when numerical breakdown occurs at step k. In this case, observe 
from Algorithm [3] that since we set tk+ik = 1, we have 

k 

Pk+1 = % - ^ UkPi 

i=\ 

with ijfc = q^{Aqk + Bpk). Therefore, we obtain 

k 

< i+Y^\Uk\\\vi\\ 

k 



k+l\ 



< 1 + ^^(11^11 + ||i?||||p,||)||p,| 



i=l 



1 + PllE 11^*11 + ll^llll^fcllEll^^ll- (29) 



j=l i=l 

Since ||s|| = 1, we can roughly estimate 

(l|Mf + ||p,+if)V^ ^ (||M||f + ||p,^if)V^ 

(l+||p,,||2)l/2 ^(l+l^)^J|p,.||2)l/2' ^ ^ 

which can be monitored during the MSOAR procedure without computing s. Note that 
\\pj\\, j = 1,2, . . . , A: + 1 are controlled and bounded in Algorithm [3l By definition (j2ip of 
Cfc, combining (j30|) with (j29l) . we conclude that c^ is bounded and furthermore is moderate 
when ||M|| is not large. 

Remark 1. Once |e^s| is considerably smaller than one, c^ -tol is a considerable overesti- 
mate on the relative residual norms in (I24p and (I25p . As a result, if Algorithm [3] breaks down 
at step k in the sense of ([20]), then ([MD and ([25]) indicate that the MSOAR and RSOAR 
methods generally find 2k eigenpairs of ([1]) at least with accuracy at the level of tol. 

Remark 2. Suppose that the convergence tolerance is ctol. Then ()20p . ()24p and ()25p 
show that if we choose toZ such that 

Cfc|e^s| • tol < ctol (31) 

then the relative residual norms of {9,y) and {6,u) definitely drop below ctol. Since Ck is 
moderate and |e^s| < 1, we propose taking 

tol = ctol (32) 

for the numerical breakdown and deflation. Numerical experiments will demonstrate that 
this choice works very well. Furthermore, they will show that such a tol is conservative and 
a bigger tol also works equally well, say, ctol = 10^^*^ and tol = 10~^. This is expected as we 
may have |e^s| < 1 considerably. 
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4.2 Implicit restarts 

As k increases, the SOAR or MSOAR method and the RSOAR method become costly and 
impractical due to storage requirement and/or computational cost. So it is generally necessary 
to restart them for practical purposes. That is, for a given maximum k, if the methods do 
not converge yet, we select a new starting vector u~^ or uf , U2 based on the information 
available to construct a better subspace Qk{^, B; u~^) or Gk{A, B; uf , Ug ) that contains richer 
information on the eigenvectors x associated with the desired eigenvalues A and compute new 
better approximate eigenpairs with respect to Qk{A, B;u~^) or Qjt{A,B;uf ,U2)- Proceed in 
such a way until the methods converge. 

For the Krylov subspace Kk{H,v) with v = [tt , 0] , if we restarted the Arnoldi process 
implicitly, an updated starting vector v'^ would have the form 

with 7 a normalizing factor such that \\v^\\ = 1 and ip{z) a filter polynomial of a certain 
degree. Obviously, for a given polynomial ip{z), v~^ does not preserve the structure of v 
any more. As a result, the new Krylov subspace ICk{H,v~^) will lose those close connections 
to the SOAR procedure presented in Section [2] and we cannot correspond K,k{H^v~^) to a 
second-order Krylov subspace. Thus, an implicit restarting scheme for the Krylov subspace 
JCk{H, v) cannot be adapted to the SOAR method and its refined version. This suggests that 
we turn to the methods based on the MSOAR procedure, so that we can make use of the 
fundamental relation (|11|) to implicitly restart the methods based on the MSOAR procedure. 

As was seen previously, the MSOAR method and the RSOAR method will find the 2k 
exact eigenpairs of QEP ([T]) if the MSOAR procedure breaks down at step k and no defiation 
occurs before step k. Recall the basic fact that the Arnoldi method and the refined Arnoldi 
method on H find the k exact eigenpairs of H if the Arnoldi process breaks down at step k. 
The implicit restarting technique updates starting vectors repeatedly to generate new Krylov 
subspaces that are expected to contain increasingly more accurate approximations to the de- 
sired m(< k) eigenvectors until the Arnoldi process approximately breaks down at step m and 
an ?7i-dimensional approximate invariant subspace of H is found. According to Theorem [3] 
and (jlip . we will have an m-dimensional approximate invariant subspace Qm-i{A,B;ui,U2) 
of QEP ([1]) when ICm{H, v) is an approximate invariant subspace of H. Fundamentally, Theo- 
rem[3]and (jlip indicate that restarting the MSOAR procedure, i.e, updating Qm{A, B; ui, U2), 
is equivalent to restarting the Arnoldi process, i.e., updating K,m{H,v). 

Under the crucial assumption that no deflation occurs. Otto [22] has proposed an implic- 
itly restarted MSOAR algorithm with exact second-order shifts suggested. There, the exact 
second-order shifts are among the 2k — m unwanted eigenvalues of the projected QEP of QEP 
([1]) onto Gk-iiA, B; ui,U2). Otto's definition means that there might be many ways of choos- 
ing the exact second-order shifts. However, he does not suggest a definite one and simply 
says that any one of the candidates can be used a shift. As will be discussed in Section [4.21 
selecting reasonable shifts among the candidates is mathematically nontrivial, and we must 
do more in order to select 'corrrect' ones. 

Given p shifts /ii, //2, • • • , jJ-p, performing p implicit shifted QR iterations on Tk yields 

(Tk - /ill) • • • (Tfe - /ipl) = VkR, 

where V^ is a A; x A; orthogonal matrix and R is upper triangular. Specifically, Vk has only p 
nonzero subdiagonals. The following results are presented in [22j . 

Theorem 6. Given p shifts //i, . . . ,/ip, perform p steps of implicit shifted QR iterations on 
Tk. Let ip{Tk) = VkRk with ipip.) = 11^=1 (/^ - f^j) o.'^d define Q^ = QkVk and r+ = V^TkVk- 
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Assume that no deflation occurs in the k{= m + p)-step MSOAR decomposition ([21 
we have an updated m-step MSOAR decomposition 
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Theorem [6] states that if there is no deflation in ()2ip then we have naturally obtained an 
m-step MSOAR procedure after p implicit shifted QR iterations are run on T^. ([33|) is then 
extended to a fc-step one from step m + 1 upwards in a standard way other than restarting 
it from scratch. 

Similar to a basic result in [24], we can prove the following theorem, which shows what 
the updated starting vector of ([33|) is. The result guides us how to suitably select shifts for 
use within implicitly restarting MSOAR and RSOAR algorithms. 



Theorem 7. R holds that 



Pt 



-^{H) 



qi 
pi 



(34) 



with ip{X) = Yl^j^ii^ — fJ-j) and r a normalizing factor. 

4.3 Selection of shifts 

In principle, we are free to use different shifts in the implicitly restarted MSOAR procedure 
described above. For linear eigenvalue problems and SVD problems, it has been shown in 
[m [13] and [T71 [18] that the better the shifts approximate some of the unwanted eigenvalues 
or singular values, the richer information on the desired eigenvectors or singular vectors 
the updated starting vector contains, so that the resulting Krylov subspaces contain more 
accurate approximations to the desired eigenvectors or singular vectors and the implicitly 
restarted algorithms are expected to converge faster. It follows from (j34p and (jlip that this 
theory works for the implicitly restarted MSOAR and RSOAR algorithms as well, that is, 
increasingly better standard Krylov subspaces lead to increasingly better modified second- 
order Krylov subspaces and vice versa. So we should choose shifts for each algorithm in the 
sense that they are best possible approximations to some of the unwanted eigenvalues of ([I|). 

We solve the projected QEP (|14p and select m Ritz values 9i as approximations to the 
desired eigenvalues. Then we may use the unwanted Ritz values as shifts, called exact second- 
order shifts in [22]. The problem is that we now have 2k Ritz values. If we used all the 
unwanted Ritz values as shifts, this would mean that we apply 2k — m>p shifts to the k x k 
matrix T^. However, for a k-step MSOAR procedure, the number of shifts in implicit restarts 
cannot exceed p. So we face the problem of how to select p < k — m reasonable shifts among 
the 2k — m candidates. 

As is known, a QEP may have two distinct eigenvalues that share the same eigenvector; 
see |27j . This means that, for QEP (J14p . if shift candidates and the m Ritz values used to 
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approximate the desired eigenvalues share the common eigenvector (s), restarting will filter 
out the information on the desired eigenvector (s). Therefore, we must abandon such shift 
candidates. Otherwise, simply selecting arbitrary p ones among 2k — m candidates, as was 
done in [22], may cause fatal consequence and make implicitly restarted algorithms fail to 
work. 

To avoid filtering out the information on the desired eigenvectors and meanwhile to 
dampen those components of the unwanted eigenvectors, we suggest a new approach to select 
reasonable shifts. We project QEP ([1]) onto the orthogonal complement of span{yi, . . . , y^} 
with respect to Qk-i{A,B;ui,U2)-, where yi,. . . ,ym are the Ritz vectors approximating the 
desired eigenvectors. Then we obtain a projected QEP of dimension p and compute its 2p 
eigenvalues. It is distinctive that these 2p eigenvalues are now approximations to some of 
the unwanted eigenvalues of QEP ([T|) since the information on xi, . . . ,Xm. has been filtered 
out from Qk-ii^, B;ui,U2)- So we can use any p ones of these 2p candidates as shifts. We 
also call them exact shifts, which are, however, different from the exact second-order shifts 
or exact shifts in Otto [22j . To be unique, we choose the p ones farthest from the Ritz values 
6i, i = 1,2, ... ,m. If we are interested in m eigenvalues nearest to a target a and/or the 
associated eigenvectors, QEP ([T|) can be equivalently transformed to a shift-invert QEP; see 
the end of Section 14.31 In this case, we select the p Ritz values among 2p candidates farthest 
from a as shifts. Such selection of shifts is motivated by the idea |171ll8j. where some of the 
shifts are taken to be unwanted Ritz values farthest from the wanted approximate singular 
values. It was argued there that this selection may better dampen those components of the 
unwanted singular vectors and amplify the components of the desired singular vectors. 

Algorithm H] computes the refined Ritz vectors ttj. Since they are generally more and 
can be much more accurate than the Ritz vectors yi [9l [161 HH] > it is possible to find better 
shifts for use within the implicitly restarted RSOAR algorithm. The first author [12^ [T3] 
has proposed certain refined shifts for the refined Arnoldi method and the refined harmonic 
Arnoldi method for linear eigenvalue problems. Later on, the same idea has been extended 
to the refined Lanczos and harmonic Lanczos bidiagonalization methods [171 [T8] . It is shown 
that the refined shifts are generally better than the corresponding exact shifts and can be 
computed efficiently and reliably. In the same spirit, we next propose certain refined shifts 
for the RSOAR method. 

Since the refined Ritz vectors Ui, i = 1,2, ... ,m are more accurate than the correspond- 
ing yi, the orthogonal complement of span{ui, . . . ,Um} with respect to Gk-ii^, B;ui,U2) 
contains richer information on the unwanted eigenvectors. As a result, the eigenvalues of 
the projected QEP of QEP ([TJ onto this complement are more accurate approximate eigen- 
values than those of the projected QEP of QEP ([1]) onto the orthogonal complement of 
span{yi, . . . ,ym} with respect to Qk-ii^i B;ui,U2). Therefore, they are generally better 
shifts than the exact shifts. We use the same approach as above to select p ones among them 
as shifts, called the refined shifts, for use within the implicitly restarted RSOAR algorithm. 

Now we show how to compute the exact and refined shifts efficiently and reliably. We 
take the refined shifts as an example. The computation of exact shifts is analogous. Recall 
Ui = QkZi, i = 1,2, . . . ,m and write Z^ = [zi, . . . , Zm]- We comment that if two columns Zi 
and Zj-|_i of Zm are complex conjugate then we replace them by their normalized real and 
imaginary parts, respectively, so that the resulting Zm is real. Make the QR decomposition 



Zm. = [Um, U±] 







where Um and U± are k x m and k x p column orthonormal matrices, respectively, and 
Rm is m X m upper triangular. We use the standard Matlab function qr.m to compute the 
decomposition in experiments. This costs 0{k^) fiops, negligible to the cost of the /s-step 
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MSOAR procedure. Then QkU± is an orthogonal basis of the orthogonal complement of 
span{ui, . . . , Um} with respect to ^fc_i(74, B; ui,U2)- It is easily seen that the projected QEP 
of the large ([T]) onto span{QfcC/_|_} is just the projected QEP of the small QEP (|18p onto 
span{U±}. Therefore, the projected QEP onto the orthogonal complement can be efficiently 
formed using 0{k^) flops. We then compute its 2p eigenvalues using 0{p^) flops and select 
p ones among them as the refined shifts. The whole cost of computing the refined shifts is 
0{k^) flops. For the exact shifts, recall the Ritz vectors yi = Qk9i, i = 1,2, ... ,m. Write 
Gm = [gi, ■ ■ ■ , dm] and replace Z^ by it. We then compute the exact shifts in the same way 
as above. 



4.4 Treatment of deflations in implicit restarts 

Previously it was assumed that no deflation occurs in implicit restarts. If deflation occurs at 
some step(s) in the MSOAR procedure, Qk must have zero column(s) and Q^ = QkVk{:, 1 : 
m) is not column orthonormal any longer and the matrix Q^ is not column orthonormal any 



longer since Qk has zero column(s), so that (j33p is not an m-step MSOAR decomposition 
any more. As a result, implicit restarting fails and is not applicable. Because of this, Otto 
|22j had to assume that deflation does not occur during the A;-step MSOAR procedure. This 
limits the applicability and generality of the algorithm. In what follows we present an effective 
approach to treat the deflation issue, so that implicit restarting works unconditionally. 



Suppose that deflation occurs at steps mi,m2. 



,mi 



< k. Then the corresponding j 



columns of Qk are zeros. Denote by Qk and Vk the matrices by deleting the zero columns 



of Qk and rows mi,m2, ■ ■ ■ , 
from which and (j2ip we get 



m. 



of Vk, respectively. Then it holds that Q~l = QkVk = QkVk, 



' A B' 




" QkVk ' 




' QkVk ' 


I 




_ PkVk _ 




_ PkVk _ 



Tu + tk+n 



Qk+i 
Pk+i 



elVk, 



(35) 



where T^ = V^TkVk- Keep in mind that Qk itself is column orthonormnal. 

Note that Vk is a {k — j) x k matrix and its row rank is k — j. For brevity, we temporarily 
assume that the first k — j columns of Vk are linearly independent, i.e., the matrix consisting 
the first k — j columns of Vk is nonsingular. Applying the Gram-Schmidt orthogonalization 
with refinement to the columns of Vk, we compute 



Vk = UkRk = [Uk-j,0] 



R 



■k-j 





Rl2 
R4 



(36) 



a variant of the QR decomposition, where Uk-j is a. {k — j) x (k — j) orthogonal matrix 
and Rk is a k X k nonsingular upper triangular matrix, and during the orthogonalization we 
set Rk{i, i) = 1, i = k — j + I, . . . ,k, that is, Rj{i, i) = 1, i = 1,2, . . . ,j. We remark that 
Uk-jRk-j is the QR decomposition of the matrix consisting of the first k — j columns of Vk. 
Then we can transform (I35D to 



A B 
I 



QkUk 
Pk^kRk 



QkUk 
PkVkRk 



RkTk ^k + ^k+lk 



Qk+1 
Pk+1 



elVkR, 



-1 



(37) 



Since Rf, is upper triangular, RkT^Rj^ is Hessenberg. Note that Vk has only p = k — m 
nonzero subdiagonals. Then the first possible nonzero entry of e^V^ is in position m and 

tk+ikelVkR',:^ = {0,...,0,^,b^) 
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with /3 = tk+ikVk{k,m)/e*^Rken 
sides of (j37|) . we obtain 



in position m. Equating the first m columns on both hand 



' A B' 

I 






= 








9m+l 
Pm+1 



(38) 



where Q^ = QkUk{'-, 1 : rn), Pm. = PkVki'--, 1 : fn)R^ with /?„ the m x m leading principal 



matrix of Rk, T^ the m x m leading principal matrix of RkT^R^ 

QkUk 



and 



Pm+1 



Jn 



with f- 



PkVkRl 



m+l'-m+lm 



+ 



Qk+1 
Pk+i 



P and /3^ 



Since Qk is column orthonormnal, Q'^ is column orthonormal too when Uk{'-, 1 : rn) is so. 
The column orthonormality of L^fc(:, 1 : m) is guaranteed whenever m < k — j, i.e., j < k — m. 
This means that Q^ is column orthonormal and no deflation occurs in (j38p provided that 
the number of deflations during the last cycle of MSOAR procedure does not exceed k — m. 
If j > /c — 771, there are m — [k — j) deflations in (jSSp and m — (k — j) zero columns among 
Q^, which correspond to the zero columns of [/(:, 1 : m). In both cases, it is easily justified 
that (Qm.)*'/m+i ~ 0- This means that ([38|) is a truly ?7i-step MSOAR procedure. So we have 
successfully cured deflations in implicit restarts and developed a robust implicit restarting 
scheme for the MSOAR procedure. 

Finally, we discuss the case that Vk has row rank k — j but its first k — j columns of Vk 
are linearly dependent. In this case, Uk in (I36p still has k — j orthonormal columns but they 
are not the first k — j ones, and the j diagonal entries ones of Rk in (j36p are not in the right 
bottom corner of Rk any more. If C/fc(:, 1 : m) has zero columns, defiation occurs in (I38p since 
Qm — QkUki'-, 1 : "rn) has zero columns, which correspond to the zero columns of C/fc(:, 1 : m). 

Having done the above, we have finally developed the following Algorithm [5l 

Algorithm 5. The implicitly restarted MSOAR type algorithms 

1. Given starting vectors qi and pi, the number m of desired eigenpairs and the integer p 
such that k = m + p, run the k-step MSOAR procedure to generate Qk- 

2. Do until convergence 

Project QEP ^ onto span{Qk} to get QEP ([M]) and solve it by the MSOAR method 
or the RSOAR method, respectively. 

Determine the convergence of m desired approximate eigenpairs {9i,yi) or m refined 
approximate eigenpairs {6i,Ui), respectively. 

3. If not converged, compute the p exact shifts or refined shifts and implicitly restart the 
MSOAR method with the p exact shifts or the RSOAR method with the p refined shifts, 
respectively. 

4. EndDo 

Algorithm [5] includes two algorithms: the implicitly restarted MSOAR algorithm with 
the exact shifts (IMSOAR) and RSOAR algorithm with the refined shifts (IRSOAR). We 
determine the convergence of a Ritz pair {9,y{= Qkd)) by the relative residual norm 

\\{e^M + 9C + K)Qkg\\ 



ll^lli + l|C||i + ll^^lli 
For the convergence of a refined Ritz pair {9,u{= Qkz)), we replace the above g by z. 
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If m eigenvalues closest to a target a are desired, we can use the shift-invert transformation 
p= xh ^it^ det{Q{a)) 7^ to transform QEP ([U to 

Q{p)x = {p^M + pC + k)x = 0, 

where M = a^M + aC + K,C = C + 2aM, K = M. Let i = -M-^C and B = -M'^K. 
In all the previous algorithms, we use such A and B. Let (p, y) be an approximate eigenpair 
(either a Ritz or refined Ritz pair) of Q{p)x = and f = Q{p)y. Then (- + a,y) is the 
corresponding approximate eigenpair of Q{X)x = {\^M + AC + K)x = 0. We easily derive 



f/p^ 



{M + C/p + k/f)y 

{a'^M + aC + K + {C + 2aM)/p + M/f)y 

{{\ + afM + (i + a)C + K)y = Q{\ + a)y 



(39) 



from which it is direct to get ||f|| from ||f|| without computing f explicitly. 

5 Numerical experiments 

In this section we report numerical examples to illustrate the performance of IMSOAR and 
IRSOAR and the superiority of IRSOAR to IMSOAR. All the experiments were run on 
an Intel(R) core(TM)2 with CPU 1.86GHz and 2GB RAM using Matlab 7.1 with e^ach = 
2.22 X 10~^^ under a Window XP system. We list CPU timings in seconds of the main parts 
in the two algorithms. The abbreviations of main parts are shown in Table [H where we also 
list the time of computing residuals of approximate eigenpairs, which may not be negligible 
when a large number of eigenpairs are required since we have to compute them directly for 
different 0's. Note that SMALL is the time of solving small QEP's plus that of computing 
refined Ritz vectors for IRSOAR. 



Table 1: Abbreviations 



TOTAL 


Total CPU time 


SOAR 


The CPU time of MSOAR procedure 


EXP 


The CPU time of exphcit projection 


IMP 


The CPU time of implicit restarting 


SMALL 


The CPU time of solving small (projected) QEP's 


RES 


The CPU time of computing residuals 


Restarts 


The number of restarts 



For each example, we used the same starting vector generated randomly in a uniform 
distribution for both algorithms. We transformed the projected QEP (J14p to the generalized 
eigenvalue problem 

-Ck —Kk 
I 





' 0g' 
g 


= e 


' Mk ■ 
/ 




' Og' 

g 



and solved it by the QZ algorithm. We can recover an eigenvector g of QEP ()14p from either 
the first k components or the last k components of [6g'^,g^]^. From the backward error 
analysis [7], it is preferable to take the first k ones if |^| > 1 and the last k ones if \0\ < 1. We 
adopted this choice. In all the examples, the breakdown or deflation tolerance is dtol and the 
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convergence tolerance is ctol. Based on (j32p and the comments followed, we took dtol = ctol 
or dtol = 100 X ctol in experiments. 

Example 1. Consider the free vibration of a string with clamped ends in a spatially 
inhomogeneous environment [29]. The equation characterizing the wave motion is described 

by 

utt + ea{x)ut = An, x S [0, vr], e > 0, 
u{t,Q) =u{t,TT) = 0. 

Approximating u{x,t) by Y12=iQkit)sin{k'Kx) and applying the Galerkin method leads to a 
second-order differential equation 

Mq{t) + Cq{t) + Kq{t) = 0, 

where q{t) = [qi{t), ...,g„(t)]^. We then get the QEP 

iX'^M + XC + K)x = 0, 

where M = (7r/2)/„, K = (7r/2)diag(j2) and 



O — \^ij ) 5 ^i 



-«J 



ea{x) sm{ix) sm{jx)dx 



We took n = 1200, e = 0.6, and a(x) = ^^(vr - xf - 201. We set a = 0.6 + 0.8i 
and were interested in 60 eigenvalues nearest to a and the corresponding eigenvectors. By 
taking dtol=10~^ , ctol=10~^^ , we tested the algorithms for various k and p and limited the 
maximum restarts to 100. Table [2] reported the results obtained by IRSOAR and IMSOAR. 







Table 2: Example 1, 


dtol=10~ 


'^, ctoh 


=10-10. 






Algorithm 


k 


P 


TOTAL 


SOAR 


EXP 


IMP 


SMALL 


RES 


Restarts 


IRSOAR 


76 


5 


195.17 


16.87 


98.22 


6.31 


62.96 


8.08 


16 


IMSOAR 


76 


5 


143.28 


17.31 


103.11 


6.73 


5.54 


8.63 


17 


IRSOAR 


76 


10 


118.98 


17.52 


57.25 


3.97 


34.05 


4.86 


9 


IMSOAR 


76 


10 


86.23 


17.38 


55.66 


3.91 


3.05 


4.77 


9 


IRSOAR 


76 


15 


107.95 


20.52 


49.77 


3.84 


28.3 


4.14 


8 


IMSOAR 


76 


15 


240.77 


46.41 


145.87 


21.69 


7.54 


11.78 


25 


IRSOAR 


86 


16 


170.55 


22.83 


75.75 


5.98 


47.55 


5.31 


9 


IMSOAR 


86 


16 


134.05 


24.84 


83.68 


6.77 


4.7 


5.61 


10 


IRSOAR 


96 


16 


165.66 


21.97 


72.06 


5.80 


55.16 


4.53 


7 


IMSOAR 


96 


16 


109.77 


21.86 


70.93 


5.84 


5.22 


4.55 


7 



We see from the table that for fixed subspace dimension k IRSOAR with larger p used 
fewer restarts and CPU timings, as expected, but it may not be the case for IMSOAR. For 
k = 76, p = 15, IMSOAR used 25 restarts, much more than those used with the smaller 
p. This is not unusual since SOAR may converge irregularly. More precisely, the residual 
norms obtained by MSOAR and its restarted versions may behave irregular even if projection 
subspaces are increasingly improved and contain more accurate approximations to the desired 
eigenvectors; see [9] for a theoretical justification. We found that from the fourteenth to the 
last restart IMSOAR with k = 76, p = 15 had one deflation at each restart and broke down 
at step 63 in the last restart, delivering an invariant subspace of dimension 62 and thus 62 
converged eigenpairs whose norms were smaller than dtol. For A; = 86 and p = 16, IRSOAR 
had one deflation at each of the last two restarts and IMSOAR had one deflation at each of 
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the last four restarts. Both algorithms broke down at the very last restart. Since a large 
number of eigenpairs are desired, computing refined Ritz vectors was quite time consuming. 
As the table indicated, explicit projection and solving the projected QEP's dominated the 
whole computational cost of IRS OAR. 

Figures [TH2] described convergence processes of two algorithms for k = 76, p = 15 and 
k = 86, p = 16, in which the left figures depicted the largest relative residual norms of 
approximate eigenpairs, while the right figures exhibited the number of deflations at each 
restart. 

It was seen that IRS OAR computed all the desired eigenvalues after eight restarts and 
there was no deflation. In contrast, IMSOAR used 25 restarts and one deflation occurred 
at each of the 14th to the 20th restarts . When we enlarged k to 86 and set p = 16, both 
algorithms converged quickly and there was no deflation during all restarts. 

From the experiments, we observed that IMSOAR was sensitive to k and p but IRSOAR 
behaved very smooth and converged faster for larger k and fixed p as well as fixed k and 
larger p. 
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Figure 1: Example 1. Left: residuals versus restarts; right: deflations versus restarts, a 
0.6 + 0.8z, k = 76, p = 15, dtol=W-^, ctol=10-^^. 
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Figure 2: Example 1. Left: residuals versus restarts; right: deflations versus restarts, a 
0.6 + 0.8z, A; = 86, p= 16, dtol=W-^, ctol=W~^^. 



We also tested two algorithms by computing six eigenvalues for dtol = ctol = 10 ^ and 
various k and p. The pairs {k,p) were (20,5), (20,8), (20,11), (25,5) and (30,5). We found 
that both algorithms converged quickly and used no more than five restarts for all pairs. We 
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found that the restarts were reduced by using the same p and larger k for both algorithms. 
There was no deflation in both algorithms for all pairs. 

Example 2. Consider the QEP arising from an ?i-degree-of- freedom damped mass-spring 
system [27]. By taking rrii = 1 and letting all the springs (respectively, dampers) have the 
same constant k (respectively, r) except ki = k^ = 2k and n = r„ = 2r, the resulting 
matrices are 



M = I, C = T-tridiag(-l,3,-l), K 



tridiag(— 1,3, — 1) 



which are very sparse. We took n = 5000, k = 5 and r = 10 and were interested in the six 
eigenvalues nearest to o" = — 13 + 0.4z and the corresponding eigenvectors. For dtol = 10~^ 
and ctol = 10"^'^, we tested IRSOAR and IMSOAR for some pairs {k,p) and limited the 
maximum restarts to 100. Table [3] lists the results, and Figures [SHU depict the convergence 
processes for two sets of parameters {k,p). 







Table 3: Example 2, dtol=10' 


-^, ctol-- 


=10-10. 






Algorithm 


k 


P 


TOTAL 


SOAR 


EXP 


IMP 


SMALL 


RES 


Restarts 


IRSOAR 


40 


15 


96.52 


24.24 


9.14 


27.17 


18.12 


7.02 


57 


IMSOAR 


40 


15 


79.66 


26.06 


10.83 


29.41 


2.73 


7.44 


60 


IRSOAR 


40 


23 


97.23 


31.75 


9.49 


26.20 


17.00 


6.44 


55 


IMSOAR 


40 


23 


91.50 


36.19 


11.04 


30.55 


2.69 


7.45 


63 


IRSOAR 


40 


31 


125.86 


46.78 


11.56 


33.81 


20.26 


8.11 


65 


IMSOAR 


40 


31 


126.52 


55.22 


13.54 


39.88 


3.52 


9.34 


78 


IRSOAR 


45 


31 


107.31 


39.61 


9.50 


28.50 


17.80 


5.80 


44 


IMSOAR 


45 


31 


104.00 


45.48 


10.84 


33.44 


3.16 


7.23 


51 


IRSOAR 


50 


31 


87.05 


26.36 


8.11 


24.75 


15.88 


4.66 


32 


IMSOAR 


50 


31 


81.41 


30.91 


9.64 


29.66 


3.05 


5.52 


38 



It can be found from Table [3] that IRSOAR used fewer restarts than IMSOAR for all 
^iven k and p. For fixed p = 31, restarts of both algorithms decreased with increasing k. 
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Figure 3: Example 2. Residuals versus restarts, dtol=W ^, ctol=W ^^. Left: k = 40, p = 15; 
right: k = 40, p = 31. 

Figures [SHU indicates that for {k = 45, p = 31) IRSOAR and IMSOAR converged irreg- 
ularly but for other {k,p) pairs the two algorithms converged quite smoothly. There was no 
deflation for both algorithms for all {k,p) pairs. 
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Figure 4: Example 2. Residuals versus restarts, dto/=10 ^,ctol=10 ^^. Left: A; = 45, p = 31; 
right: /c = 50, p = 31. 



Example 3. This problem comes from [2]. It is a nonlinear eigenvalue problem modeling 
a radio-frequency gun cavity that is of the form 

r(A).T = [K - \M + i(A - alY^^Wi + i{X - alfl'^W'^x = 0, 

where M, i^, VFi, W2 are real symmetric matrices of size 9956 x 9956. From these matrices, 
we constructed the following QEP of the form 

{X^W2 + XM + K)x = 0. 

We used IRSOAR and IMSOAR for dtol = 10"^, ctol = 10"^ to compute the 24 eigen- 
values nearest to o" = 0.5 -|- 0.5i and the associated eigenvectors. We stopped the algorithms 
after 50 restarts were used. Table H] reported the results. 



Table 4: Example 3, dtol=lQ-^ , cto/=10"^. 



Algorithm 


k 


P 


TOTAL 


SOAR 


EXP 


IMP 


SMALL 


RES 


Restarts 


IRSOAR 


40 


5 


22.41 


15.93 


1.4 


1.05 


1.36 


2.05 


1 


IMSOAR 


40 


5 


- 


- 


- 


- 


- 


- 


50 


IRSOAR 


40 


10 


21.39 


15.27 


1.49 


0.95 


1.43 


1.73 


1 


IMSOAR 


40 


10 


- 


- 


- 


- 


- 


- 


50 


IRSOAR 


40 


15 


22.33 


16.34 


1.51 


0.922 


1.15 


1.91 


1 


IMSOAR 


40 


15 


- 


- 


- 


- 


- 


- 


50 


IRSOAR 


50 


10 


19.05 


15.47 


1.07 





1.00 


1.17 





IMSOAR 


50 


10 


- 


- 


- 


- 


- 


- 


50 



It is remarkable to see that for /c = 40 and p = 5, 10 and 15, IRSOAR used only one restart 
while IMSOAR failed to converge after 50 restarts. So the IRSOAR is much more efficient 
than IMSOAR. We have also observed that there was no deflation for both algorithms. If 
we instead took a bigger A; = 50 and p = 10, we found that IRSOAR computed the desired 
eigenpairs without restart but IMSOAR still did not converge. So for this problem, IRSOAR 
worked very well but IMSOAR completely failed for given {k,p). 

To better look into the convergence behavior. Figure [5] depicted the convergence curves 
of both algorithms for A; = 40, p = 10 and A; = 40, p = 15. We see that for IMSOAR the 
maximum relative residuals of Ritz pairs stayed roughly between 10~^ ~ 10~^ and did not 
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decrease further as restarts proceeded. Furthermore, keep in mind that in the first cycle 
IRSOAR and IMSOAR used the same starting vector and computed approximate eigenpairs 
with respect to the same subspace. It is strikingly seen from Figure \5\ that in the first cycle 
the maximum residual norms obtained by RSOAR improved those obtained by MSOAR by 
at least three orders, a very big gain. Experimentally, this confirms the theory of |9j and 
indicates that one can benefit much from the refined projection principle. 
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Figure 5: Example 3. Residuals versus restarts. Left: k = 40, p = 10; right: k = 40, p = 15. 



Example 4. In this example, the QEP Q{X)x = (A^Af + AC + (1 + ifi)K)x = arises 
in a model of a concrete structure supporting a machine assembly [2l [5] . The matrices have 
dimension 2472, where M is real diagonal and of low rank, and C is the viscous damping 
matrix and is pure imaginary and diagonal. The factor 1 + i/i adds uniform hysteretic 
damping. The default is // = 0.04. 

We ran IRSOAM and IMSOAR to compute the 30 eigenvalues of Q{X)x = nearest to 
0" = 34 + lOi and the corresponding eigenvectors. Table [5] reported the results. 







Table 5: Example 4, dtol=10 


"^, ctol 


=10-8. 






Algorithm 


k 


P 


TOTAL 


SOAR 


EXP 


IMP 


SMALL 


RES 


Restarts 


IRSOAR 


42 


5 


19.13 


6.48 


1.16 


2.25 


4.31 


2.48 


9 


IMSOAR 


42 


5 


15.36 


6.48 


1.18 


2.28 


0.64 


3.02 


9 


IRSOAR 


42 


8 


13.97 


6.26 


0.72 


1.34 


2.79 


1.13 


5 


IMSOAR 


42 


8 


11.30 


5.87 


0.64 


1.23 


0.38 


1.70 


5 


IRSOAR 


42 


11 


10.94 


5.88 


0.54 


1.02 


1.45 


1.75 


4 


IMSOAR 


42 


11 


9.38 


5.90 


0.55 


0.97 


0.31 


1.39 


4 


IRSOAR 


45 


8 


11.77 


5.62 


0.68 


1.17 


2.22 


1.23 


4 


IMSOAR 


45 


8 


9.75 


5.65 


0.64 


1.11 


0.39 


1.59 


4 


IRSOAR 


50 


8 


11.78 


5.41 


0.66 


1.05 


2.48 


1.50 


3 


IMSOAR 


50 


8 


9.22 


5.47 


0.64 


1.03 


0.44 


1.42 


3 



It is seen that this problem is easy to solve by both IRSOAR and IMSOAR. The algorithms 
converged smoothly and quickly. For fixed p = 8 and k = 42, 45, 50, both algorithms used 
fewer restarts with increasing k. These phenomena are in accordance with the theory since 
subspaces of increasing dimensions k + p should generally contain more information on the 
desired eigenvectors. For p = 5, IRSOAR and IMSOAR had one defiation at the 39th, 40th 
and 41th restarts, and they broke down at the 41th restart. For p = 8, both algorithms had 
one deflation at each of the last four restarts, and broke down at the 42th restart. 
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Figure 6: Example 4. Left: Residuals versus restarts; right: deflations versus restarts. 
k = A2,p = 8. 



Figure [6] depicted the convergence processes and deflation details with k = 42, p = 8, 
from which we see that both IRSOAR and IMSOAR had four deflations at the last restart. In 
addition, the figure illustrates that although IRSOAR and IMSOAR both used five restarts 
to achieve the convergence, IRSOAR gave smaller residuals, i.e., more accurate approximate 
eigenpairs, at each restart. 

Example 5. We consider the damped vibration mode of an acoustic fluid confined in 
a cavity with absorbing walls capable of dissipating acoustic energy [9]. We take the same 
geometrical data as in [9] and more properties on this problem can be found there. The QEP 
is 

X^MuU + (a + \/3)Au + Kuu = 0, 

where a = 5 x 10^ N/m^,/3 = 200Ns/m^. The dimension of this QEP is n = 46548. 

We ran IRSOAM and IMSOAR to compute the 8 eigenvalues nearest to cj = 25 + 18i and 
the corresponding eigenvectors. We took dtol = 10~^, ctol = 10"^*^ in experiments. Table [6] 
reported the results. 
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1.77 


1.31 
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35 


20 


30.15 


26.66 


1.46 





0.05 


1.34 






The table clearly indicates that IRSOAR solved the problem very successfully and, re- 
garding restarts, was several times (from nearly one to 17) faster than IMSOAR except for 
k = 35, in which case both algorithms found the desired eigenpairs without restart. We fur- 
ther observe that for fixed k, IRSOAR used fewer restarts with increasing p while IMSOAR 
converged quite irregularly with varying p: If p is increased from 10 to 15, the number of 
restarts was reduced from 9 to 4; if we changed p to 20, the number of restarts was increased 
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to 18. This may be due to the intrinsically irregular convergence behavior of the Rayleigh- 
Ritz method for QEP [9] It is also seen that for fixed p = 20 if we changed k from 30 to 
35 then IMSOAR drastically converged without restart. These results demonstrate that, un- 
like IRSOAR, IMSOAR can be quite sensitive to parameter choices, so that the convergence 
behavior of IRSOAR is more pronounced than that of IRSOAR. 
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Figure 7: Example 5. Left: Residuals versus restarts; right: deflations versus restarts. 
k = 30,p = 10. 
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Figure 8: Example 5. Left: Residuals versus restarts; right: deflations versus restarts. 
k = 30,p = 20. 



Figures [THE] described the convergence processes of two algorithms for k = 30, p = 10 and 
k = 30, p = 20. Obviously, IRSOAR converged smoothly and quickly and there was no defla- 
tion occurred. In contrast, IMSOAR converged irregularly for all given pairs. Furthermore, 
we observed that for A; = 30 and p = 10, IMSOAR had three deflations at steps 23, 25 and 
26 and it broke down at the 26th step in the last restart; for k = 30, p = 20, IMSOAR had 
four deflations that occurred at steps 12,13,14 and 15 and broke down at the 15th step at 
the last restart, delivering 15 converged eigenpairs. 

6 Conclusion 

Based on the refined projection principle and the SOAR procedure, we have proposed a refined 
second-order Arnoldi (RSOAR) method that has the same structure preserving properties as 
the SOAR method. We have shown that the original SOAR method cannot be restarted 
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effectively and instead the MSOAR procedure is needed. We have estabhshed upper bounds 
for residual norms obtained by the MSOAR and RSOAR algorithms. From them, we have 
proposed reasonable criteria on numerical deflation and breakdown. The criteria ensure 
that the algorithms converge within the prescribed accuracy and meanwhile maintain the 
numerical stability of MSOAR procedure. For practical purposes, we have extended the 
implicit restarting technique to the MSOAR and RSOAR methods, developed their implicitly 
restarted algorithms IRSOAR and IMSOAR and proposed best possible shifts for each of 
them. We have shown how to compute the shifts efficiently and reliably. 

Our final contribution is that, in order to overcome possible deflation and make implicit 
restarting unconditionally applicable, we have proposed an effective approach to handle the 
deflation issue. Numerical examples have illustrated that IRSOAR and IMSOAR generally 
works well but the former is more robust and can be considerably more efficient than the 
latter. 
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